# This file and accompanying functions are renamed with domestic RC added to 
# a CF on BLP data
#
rm(list=ls())
library(HeadR)
library(SQUAREM)
library(fixest)
library(tictoc)
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_domestic_functions.R")
source("Rfiles/CF_MMX_functions.R")
source("Rfiles/CF_MMX_MLOG_functions.R")
library(Rcpp)
# create the Rcpp function we need for computing equilibrium
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp")
#parallel processing set up for doParallel, doRNG
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
#
n.hh <- 100000 # careful, this is slow, but works
U0 <- 1  # needed in prob.mat(), but not prob.mat.namedxvars()
n.reps <- 1 # for richsubs (no need to do a lot if a lot of consumers.)
sigmadomesticpct <- 0 # either  0, 72 or 172
#
# This is intended to choose which country is targeted
source(file = "Rfiles/CF_BLPdata_brandiso.R")
target.list <- brands[iso!="USA",unique(iso)] # all non US brands
rm(brands)
#
#
# Here choose which set of parameters to use. 
param.all <- fread("Data/BLP99/pyblp_domestic_param.csv") # estimated using C&G package
param.list <- names(param.all[,par_abbrev:=NULL])
param.list 
#add pub_param (the BLP 99 parameters)
param.list <- c("pub_param","Domestic_pub_param",param.list)
#
#
#Loop over param list, conducting counterfactuals ====
tic()
for(param.chosen in param.list ) {
# This is used in the counterfactuals : number of characteristics (outside the constant)
if(param.chosen %like% "^Domestic"){
n.x <- 5 #THIS NEEDS TO BE 5 WHEN ADDING DOMESTIC, 4 otherwise 
} else n.x <- 4
#
###
# What is the average own elasticity? ====
###
setwd("~/Dropbox/BLP_new")
source("Rfiles/CF_BLPdata_domestic_getdata.R") # get original data: original published_param.csv is 1995 and should not be changed (used in other files)
#
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] # super assignment needed here because of scope for DGP.BLP
BLP.meanownelas
#
###
# Implement CF on BLP ====
###
#
source("Rfiles/CF_BLPdata_domestic_getdata.R") # get original data
year.gph <- 1990
tar <- 0 # original tariff (to be raised)
tar.increase <- 0.1 # the actual increase in tariff rate
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix 
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#do a settings file only for the speed 
settings.var <- data.frame(speed =c(0.3,0.3,0.3))
#
# Run the calibration to get xi and mc
year.monte <- year.gph
objective <- "CF"
market <- 1
#
#Run sims for the 3 versions
#EHA:
set.seed(seedn)
CF_BLPdata.simulv1 <- CF_BLPdata.simul(vsim=1,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv2 <- CF_BLPdata.simul(vsim=2,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.simul(vsim=3,cf.meth="EHA")
#
CF_BLPdata.simulallv <- rbind(CF_BLPdata.simulv1,CF_BLPdata.simulv2,CF_BLPdata.simulv3) 
saveRDS(CF_BLPdata.simulallv,paste0("Tables/CF_BLPdata/CF_BLPdata_simul_",param.chosen,".rds"))
}
toc()

#now collect the RDS of output from above and create tex files
#Make the tables in tex ====
for(param.chosen in param.list ) {
CF_BLPdata.simulallv <- readRDS(paste0("Tables/CF_BLPdata/CF_BLPdata_simul_",param.chosen,".rds"))
CF_BLPdata.simulallv[, v_name:="BLP"]
CF_BLPdata.simulallv[v==1, v_name:="Logit"]
CF_BLPdata.simulallv[v==2, v_name:=" $\\beta $ het."]
CF_BLPdata.simulallv[v_name=="BLP", v_table:=1]
CF_BLPdata.simulallv[v_name=="Logit", v_table:=2]
CF_BLPdata.simulallv[v_name==" $\\beta $ het.",  v_table:=3]
#select and order the output for the tables
CF_BLPdata.table <- CF_BLPdata.simulallv[,list(v_table,v_name,meanownelas=-meanownelas,
                                               meansuperelas,pt,pte,
                                               chg_true_qty,chg_ces_qty)]
#
#save a latex result file
setorder(CF_BLPdata.table, v_table)
CF_BLPdata.table[,c("v_table"):=NULL]
#
CF_BLPdata.table[,output := texout(.SD,digits =2)]
if(length(target.list) == 1){
  writeLines(CF_BLPdata.table$output,paste0("Tables/CF_BLPdata/target_",target.list,"/CF_BLPdata_results_",param.chosen,sigmadomesticpct,".tex"))
} else writeLines(CF_BLPdata.table$output,paste0("Tables/CF_BLPdata/CF_BLPdata_results_",param.chosen,sigmadomesticpct,".tex"))
}
